
##########################################################################################################
##      Wen, Fangqi, Erik H. Wang, and Michael Hout. 2024.                                              ##
##      "Social Mobility in the Tang Dynasty as the Imperial Examination Rose and                       ##
##      Aristocratic Family Pedigree Declined, 618-907 CE."                                             ##
##      Proceedings of the National Academy of Sciences 121(4): e2305564121.                            ##
##      Replication Codes                                                                               ##
##      Jan 18, 2024                                                                                    ##
##########################################################################################################


rm(list=ls())
library(ggpubr)
library(ggplot2)
library(fixest)
library(gridExtra)

## please set your directory right

load("PNAS_Tang_Mobility_Replication.RData")


### Figure 1. The Rubbing of a Tomb Epitaph of a Deceased Elite

### Figure 2. The Distribution of Office Rank Among the Tang Elites

########
People$Rank10 <- ifelse(People$Rank_New<=19 & People$Rank_New>=18, 9,
                        ifelse(People$Rank_New<=17 & People$Rank_New>=16, 8,
                        ifelse(People$Rank_New<=15 & People$Rank_New>=14, 7,
                        ifelse(People$Rank_New<=13.5 & People$Rank_New>=12, 6,
                        ifelse(People$Rank_New<=11.5 & People$Rank_New>=10, 5,
                        ifelse(People$Rank_New<=9.5 & People$Rank_New>=8, 4,
                        ifelse(People$Rank_New<=7.5 & People$Rank_New>=6, 3,
                        ifelse(People$Rank_New<=5.5 & People$Rank_New>=4, 2,
                        ifelse(People$Rank_New<=3.5 & People$Rank_New>=2, 1,
                        ifelse(People$Rank_New<=1 & People$Rank_New>=0, 0, NA))))))))))
summary(as.factor(People$Rank10))

People$SubRank <- ifelse(People$Rank_New==19 | People$Rank_New==17 | People$Rank_New==15,
                         1, 
                         ifelse(People$Rank_New==18 | People$Rank_New==16 | People$Rank_New==14,
                         2,
                         ifelse(People$Rank_New==13.5 | People$Rank_New==11.5,
                         3,
                         ifelse(People$Rank_New==13 | People$Rank_New==11,
                         4,
                         ifelse(People$Rank_New==12.5 | People$Rank_New==10.5,
                         5,
                         ifelse(People$Rank_New==12 | People$Rank_New==10,
                         6,
                         ifelse(People$Rank_New==9.5 | People$Rank_New==7.5,
                         7,
                         ifelse(People$Rank_New==9 | People$Rank_New==7,
                         8,
                         ifelse(People$Rank_New==8.5 | People$Rank_New==6.5,
                         9,
                         ifelse(People$Rank_New==8 | People$Rank_New==6,
                         10,
                         ifelse(People$Rank_New==5.5 | People$Rank_New==3.5,
                         11,
                         ifelse(People$Rank_New==5 | People$Rank_New==3,
                         12,
                         ifelse(People$Rank_New==4.5 | People$Rank_New==2.5,
                         13,
                         ifelse(People$Rank_New==4 | People$Rank_New==2,
                         14,
                         ifelse(People$Rank_New==1, 15, 16)))))))))))))))

summary(People$SubRank)



### Save legends for all the subrank categories

get_legend <- function(p) {
  g <- ggplotGrob(p)
  k <- which(g$layout$name=="guide-box")
  g$grobs[[k]]
}


figure2 <- ggplot(data=subset(People, People$Rank10==0), aes(x=as.factor(Rank10), 
                              fill=as.factor(SubRank), color=as.factor(SubRank))) + 
  geom_bar(aes(y = (..count..)/sum(..count..))) +
  theme_bw() + 
  xlab("Office Rank") + ylab("Proportion") +
  scale_x_discrete(labels=c("0"="No Office\n or No Rank", "1"="9", "2"="8",
                            "3"="7", "4"="6", "5"="5", "6"="4", "7"="3", "8"="2", "9"="1")) +
  scale_fill_manual(name = "No Rank / Office",        
                    labels=c("15"="No Rank", "16"="No Office"),
                    values=c("#252525", "#969696")) +
  scale_color_manual(name = "No Rank / Office", 
                     labels=c("15"="No Rank", "16"="No Office"),
                     values=c("#252525", "#969696")) +
  theme(axis.title.x = element_text(face = "bold", size=20),
        axis.title.y = element_text(face = "bold", size=20), 
        plot.margin=unit(c(0.5, 0.5, 0.5, 0.5),"cm"),
        axis.text.x=element_text(face = "bold", size=20),
        axis.text.y=element_text(face = "bold", size=20),
        legend.title=element_text(size=16, face = "bold"),
        legend.text=element_text(size=16)) 
figure2

no_rank_no_office <- get_legend(figure2)



figure2 <- ggplot(data=subset(People, People$Rank10==1), aes(x=as.factor(Rank10), 
                              fill=as.factor(SubRank), color=as.factor(SubRank))) + 
  geom_bar(aes(y = (..count..)/sum(..count..))) +
  theme_bw() + 
  xlab("Office Rank") + ylab("Proportion") +
  scale_x_discrete(labels=c("0"="No Office\n or No Rank", "1"="9", "2"="8",
                            "3"="7", "4"="6", "5"="5", "6"="4", "7"="3", "8"="2", "9"="1")) +
  scale_fill_manual(name = "Low Rank", 
                    labels=c("11"="Principal Upper", "12"="Principal Lower",
                             "13"="Junior Upper", "14"="Junior Lower"),
                    values=c("#08519c", "#3182bd", "#6baed6", "#bdd7e7")) +
  scale_color_manual(name = "Low Rank", 
                     labels=c("11"="Principal Upper", "12"="Principal Lower",
                              "13"="Junior Upper", "14"="Junior Lower"),
                     values=c("#08519c", "#3182bd", "#6baed6", "#bdd7e7")) +
  theme(axis.title.x = element_text(face = "bold", size=20),
        axis.title.y = element_text(face = "bold", size=20), 
        plot.margin=unit(c(0.5, 0.5, 0.5, 0.5),"cm"),
        axis.text.x=element_text(face = "bold", size=20),
        axis.text.y=element_text(face = "bold", size=20),
        legend.title=element_text(size=16, face = "bold"),
        legend.text=element_text(size=16)) 
figure2

low_rank <- get_legend(figure2)



figure2 <- ggplot(data=subset(People, People$Rank10==3), aes(x=as.factor(Rank10), 
                             fill=as.factor(SubRank), color=as.factor(SubRank))) + 
  geom_bar(aes(y = (..count..)/sum(..count..))) +
  theme_bw() + 
  xlab("Office Rank") + ylab("Proportion") +
  scale_x_discrete(labels=c("0"="No Office\n or No Rank", "1"="9", "2"="8",
                            "3"="7", "4"="6", "5"="5", "6"="4", "7"="3", "8"="2", "9"="1")) +
  scale_fill_manual(name = "Middle Rank", 
                    labels=c("7"="Principal Upper", "8"="Principal Lower",
                             "9"="Junior Upper", "10"="Junior Lower"),
                    values=c("#006d2c", "#31a354", "#74c476", "#bae4b3")) +
  scale_color_manual(name = "Middle Rank", 
                     labels=c("7"="Principal Upper", "8"="Principal Lower",
                              "9"="Junior Upper", "10"="Junior Lower"),
                     values=c("#006d2c", "#31a354", "#74c476", "#bae4b3")) +
  theme(axis.title.x = element_text(face = "bold", size=20),
        axis.title.y = element_text(face = "bold", size=20), 
        plot.margin=unit(c(0.5, 0.5, 0.5, 0.5),"cm"),
        axis.text.x=element_text(face = "bold", size=20),
        axis.text.y=element_text(face = "bold", size=20),
        legend.title=element_text(size=16, face = "bold"),
        legend.text=element_text(size=16)) 
figure2

middle_rank <- get_legend(figure2)



figure2 <- ggplot(data=subset(People, People$Rank10==5), aes(x=as.factor(Rank10), 
                              fill=as.factor(SubRank), color=as.factor(SubRank))) + 
  geom_bar(aes(y = (..count..)/sum(..count..))) +
  theme_bw() + 
  xlab("Office Rank") + ylab("Proportion") +
  scale_x_discrete(labels=c("0"="No Office\n or No Rank", "1"="9", "2"="8",
                            "3"="7", "4"="6", "5"="5", "6"="4", "7"="3", "8"="2", "9"="1")) +
  scale_fill_manual(name = "High Rank", 
                    labels=c("3"="Principal Upper", "4"="Principal Lower",
                             "5"="Junior Upper", "6"="Junior Lower"),
                    values=c("#e6550d", "#fd8d3c", "#fdbe85", "#fdd0a2")) +
  scale_color_manual(name = "High Rank", 
                     labels=c("3"="Principal Upper", "4"="Principal Lower",
                              "5"="Junior Upper", "6"="Junior Lower"),
                     values=c("#e6550d", "#fd8d3c", "#fdbe85", "#fdd0a2")) +
  theme(axis.title.x = element_text(face = "bold", size=20),
        axis.title.y = element_text(face = "bold", size=20), 
        plot.margin=unit(c(0.5, 0.5, 0.5, 0.5),"cm"),
        axis.text.x=element_text(face = "bold", size=20),
        axis.text.y=element_text(face = "bold", size=20),
        legend.title=element_text(size=16, face = "bold"),
        legend.text=element_text(size=16)) 
figure2

high_rank <- get_legend(figure2)


figure2 <- ggplot(data=subset(People, People$Rank10==9), aes(x=as.factor(Rank10), 
                              fill=as.factor(SubRank), color=as.factor(SubRank))) + 
  geom_bar(aes(y = (..count..)/sum(..count..))) +
  theme_bw() + 
  xlab("Office Rank") + ylab("Proportion") +
  scale_x_discrete(labels=c("0"="No Office\n or No Rank", "1"="9", "2"="8",
                            "3"="7", "4"="6", "5"="5", "6"="4", "7"="3", "8"="2", "9"="1")) +
  scale_fill_manual(name = "Elite Rank", 
                    labels=c("1"="Principal", "2"="Junior"),
                    values=c("#ae017e", "#d7b5d8")) +
  scale_color_manual(name = "Elite Rank", 
                     labels=c("1"="Principal", "2"="Junior"),
                     values=c("#ae017e", "#d7b5d8")) +
  theme(axis.title.x = element_text(face = "bold", size=20),
        axis.title.y = element_text(face = "bold", size=20), 
        plot.margin=unit(c(0.5, 0.5, 0.5, 0.5),"cm"),
        axis.text.x=element_text(face = "bold", size=20),
        axis.text.y=element_text(face = "bold", size=20),
        legend.title=element_text(size=16, face = "bold"),
        legend.text=element_text(size=16)) 
figure2

elite_rank <- get_legend(figure2)



p <- ggplot(data=subset(People, !is.na(People$Rank_New)), aes(x=as.factor(Rank10), 
                        fill=as.factor(SubRank), color=as.factor(SubRank))) + 
  geom_bar(aes(y = (..count..)/sum(..count..))) +
  theme_bw() + 
  xlab("Office Rank") + ylab("Proportion") +
  ylim(0, 0.32) +
  scale_x_discrete(labels=c("0"="No Office\n or No Rank", "1"="9", "2"="8",
                            "3"="7", "4"="6", "5"="5", "6"="4", "7"="3", "8"="2", "9"="1")) +
  scale_fill_manual(name = "", 
                    labels=c("1"="Elite Rank: Principal", "2"="Elite Rank: Junior", 
                             "3"="High Rank: Principal Upper", "4"="High Rank: Principal Lower",
                             "5"="High Rank: Junior Upper", "6"="High Rank: Junior Lower",
                             "7"="Middle Rank: Principal Upper", "8"="Middle Rank: Principal Lower",
                             "9"="Middle Rank: Junior Upper", "10"="Middle Rank: Junior Lower",
                             "11"="Low Rank: Principal Upper", "12"="Low Rank: Principal Lower",
                             "13"="Low Rank: Junior Upper", "14"="Low Rank: Junior Lower",
                             "15"="No Rank", "16"="No Office"),
                    values=c("#ae017e", "#d7b5d8", "#e6550d", "#fd8d3c", "#fdbe85", "#fdd0a2",
                             "#006d2c", "#31a354", "#74c476", "#bae4b3",
                             "#08519c", "#3182bd", "#6baed6", "#bdd7e7",
                             "#252525", "#969696")) +
  scale_color_manual(name = "", 
                     labels=c("1"="Elite Rank: Principal", "2"="Elite Rank: Junior", 
                              "3"="High Rank: Principal Upper", "4"="High Rank: Principal Lower",
                              "5"="High Rank: Junior Upper", "6"="High Rank: Junior Lower",
                              "7"="Middle Rank: Principal Upper", "8"="Middle Rank: Principal Lower",
                              "9"="Middle Rank: Junior Upper", "10"="Middle Rank: Junior Lower",
                              "11"="Low Rank: Principal Upper", "12"="Low Rank: Principal Lower",
                              "13"="Low Rank: Junior Upper", "14"="Low Rank: Junior Lower",
                              "15"="No Rank", "16"="No Office"),
                     values=c("#ae017e", "#d7b5d8", "#e6550d", "#fd8d3c", "#fdbe85", "#fdd0a2",
                              "#006d2c", "#31a354", "#74c476", "#bae4b3",
                              "#08519c", "#3182bd", "#6baed6", "#bdd7e7",
                              "#252525", "#969696")) +
  theme(axis.title.x = element_text(face = "bold", size=20),
        axis.title.y = element_text(face = "bold", size=20), # angle = 90, 
        plot.margin=unit(c(0.5, 0.5, 0.5, 0.5),"cm"),
        axis.text.x=element_text(face = "bold", size=20),
        axis.text.y=element_text(face = "bold", size=20),
        legend.position = "none",
        legend.title=element_text(size=16, face = "bold"),
        legend.text=element_text(size=16)) 
p



elite_rank <- as_ggplot(elite_rank)
high_rank <- as_ggplot(high_rank)
middle_rank <- as_ggplot(middle_rank)
low_rank <- as_ggplot(low_rank)
no_rank_no_office <- as_ggplot(no_rank_no_office)


legend <- grid.arrange(elite_rank, high_rank, middle_rank, low_rank, no_rank_no_office, ncol=1)


### Figure 2 with legends
grid.arrange(arrangeGrob(p, legend, widths = c(0.8, 0.2)))
# (You might need to adjust the space between some rank categories.)



### Table 1. Descriptive statistics

## Before and After Wu Zetian's Rule
Pre_Empress <- subset(People, People$After_Empress==0)
Post_Empress <- subset(People, People$After_Empress==1)


full <- People[, c("branch", "FRank", "Exam", "Rank_New", "GPRank", "MarriageFull", 
                   "Bieji", "Modern.province", "Pre.TangRegion_alt", "decade")]
full <- na.omit(full)
full <- full[, 1:4]

Pre <- Pre_Empress[, c("branch", "FRank", "Exam", "Rank_New", "GPRank", "MarriageFull", 
                       "Bieji", "Modern.province", "Pre.TangRegion_alt", "decade")]
Pre <- na.omit(Pre)
Pre <- Pre[, 1:4]

Post <- Post_Empress[, c("branch", "FRank", "Exam", "Rank_New", "GPRank", "MarriageFull", 
                         "Bieji", "Modern.province", "Pre.TangRegion_alt", "decade")]
Post <- na.omit(Post)
Post <- Post[, 1:4]


descriptive_mean <- as.data.frame(apply(full, 2, function (x) mean(x, na.rm = T)))
descriptive_sd <- as.data.frame(apply(full, 2, function (x) sd(x, na.rm = T)))
descriptive_min <- as.data.frame(apply(full, 2, function (x) min(x, na.rm = T)))
descriptive_max <- as.data.frame(apply(full, 2, function (x) max(x, na.rm = T)))
descriptive_n <- as.data.frame(apply(full, 2, function (x) length(x[!is.na(x)])))
descriptive_full <- cbind.data.frame(round(descriptive_mean, digits=3), 
                                     round(descriptive_sd, digits=3), 
                                     round(descriptive_min, digits=3), 
                                     round(descriptive_max, digits=3),
                                     descriptive_n)

colnames(descriptive_full) <- c("Mean", "SD", "Min", "Max", "N")

##########
descriptive_mean <- as.data.frame(apply(Pre, 2, function (x) mean(x, na.rm = T)))
descriptive_sd <- as.data.frame(apply(Pre, 2, function (x) sd(x, na.rm = T)))
descriptive_min <- as.data.frame(apply(Pre, 2, function (x) min(x, na.rm = T)))
descriptive_max <- as.data.frame(apply(Pre, 2, function (x) max(x, na.rm = T)))
descriptive_n <- as.data.frame(apply(Pre, 2, function (x) length(x[!is.na(x)])))
descriptive_Pre <- cbind.data.frame(round(descriptive_mean, digits=3), 
                                    round(descriptive_sd, digits=3), 
                                    round(descriptive_min, digits=3), 
                                    round(descriptive_max, digits=3),
                                    descriptive_n)

colnames(descriptive_Pre) <- c("Mean", "SD", "Min", "Max", "N")

##########
descriptive_mean <- as.data.frame(apply(Post, 2, function (x) mean(x, na.rm = T)))
descriptive_sd <- as.data.frame(apply(Post, 2, function (x) sd(x, na.rm = T)))
descriptive_min <- as.data.frame(apply(Post, 2, function (x) min(x, na.rm = T)))
descriptive_max <- as.data.frame(apply(Post, 2, function (x) max(x, na.rm = T)))
descriptive_n <- as.data.frame(apply(Post, 2, function (x) length(x[!is.na(x)])))
descriptive_Post <- cbind.data.frame(round(descriptive_mean, digits=3), 
                                     round(descriptive_sd, digits=3), 
                                     round(descriptive_min, digits=3), 
                                     round(descriptive_max, digits=3),
                                     descriptive_n)

colnames(descriptive_Post) <- c("Mean", "SD", "Min", "Max", "N")


descriptive_table <- cbind.data.frame(descriptive_full, descriptive_Pre, descriptive_Post)
descriptive_table



### Table 2 ###

full <- People[, c("Rank_New", "branch", "FRank", "Exam", "GPRank", "MarriageFull", "Bieji", 
                   "Modern.province", "Pre.TangRegion_alt", "decade")]
Pre <-  Pre_Empress[, c("Rank_New", "branch", "FRank", "Exam", "GPRank", "MarriageFull", "Bieji", 
                        "Modern.province", "Pre.TangRegion_alt", "decade")]
Post <- Post_Empress[, c("Rank_New", "branch", "FRank", "Exam", "GPRank", "MarriageFull", "Bieji", 
                         "Modern.province", "Pre.TangRegion_alt", "decade")]


### Table 2 Panel A

m1 <- feols(Rank_New ~ branch + Exam + FRank + GPRank + MarriageFull +
              Bieji  | Modern.province^decade +
              Pre.TangRegion_alt^decade, 
            data = Pre)
m1_clse <- summary(m1, cluster = c("decade"))
summary(m1_clse)

m2 <- feols(Rank_New ~ branch + Exam + FRank + GPRank + MarriageFull +
              Bieji  | Modern.province^decade +
              Pre.TangRegion_alt^decade, 
            data = Post)
m2_clse <- summary(m2, cluster = c("decade"))
summary(m2_clse)

m3 <- feols(Rank_New ~ branch + Exam + FRank + GPRank + MarriageFull +
              Bieji  | Modern.province^decade +
              Pre.TangRegion_alt^decade, 
            data = People)
m3_clse <- summary(m3, cluster = c("decade"))
summary(m3_clse)

m4 <- feols(Rank_New ~ branch*Exam + FRank*Exam + GPRank + MarriageFull +
              Bieji  | Modern.province^decade +
              Pre.TangRegion_alt^decade, 
            data = People)
m4_clse <- summary(m4, cluster = c("decade"))
summary(m4_clse)

etable(m1_clse, m2_clse, m3_clse, m4_clse, digits = 3,
       keep=c("branch", "Exam", "FRank"),
       signif.code=c("***"=0.001, "**"=0.01, "*"=0.05))

### Table 2 Panel B

m5 <- feols(Exam ~ branch + FRank + GPRank + MarriageFull +
              Bieji  | Modern.province^decade +
              Pre.TangRegion_alt^decade, 
            data = subset(Pre, !is.na(Pre$Rank_New)))
m5_clse <- summary(m5, cluster = c("decade"))
summary(m5_clse)

m6 <- feols(Exam ~ branch + FRank + GPRank + MarriageFull +
              Bieji  | Modern.province^decade +
              Pre.TangRegion_alt^decade, 
            data = subset(Post, !is.na(Post$Rank_New)))
m6_clse <- summary(m6, cluster = c("decade"))
summary(m6_clse)

m7 <- feols(Exam ~ branch + FRank + GPRank + MarriageFull +
              Bieji  | Modern.province^decade +
              Pre.TangRegion_alt^decade, 
            data = subset(People, !is.na(People$Rank_New)))
m7_clse <- summary(m7, cluster = c("decade"))
summary(m7_clse)

etable(m5_clse, m6_clse, m7_clse, digits = 3,
       keep=c("branch", "Exam", "FRank"),
       signif.code=c("***"=0.001, "**"=0.01, "*"=0.05))



### Figure 3 ###
plot_model <- feols(Rank_New ~ branch + 
                      branch:decade +
                      branch:I(decade^2) + 
                      Exam + 
                      Exam:decade +
                      Exam:I(decade^2) +
                      FRank +
                      FRank:decade + 
                      FRank:I(decade^2) +
                      GPRank +
                      GPRank:decade + 
                      GPRank:I(decade^2) +
                      MarriageFull + 
                      MarriageFull:decade +
                      MarriageFull:I(decade^2) +
                      Bieji | Modern.province^decade +
                      Pre.TangRegion_alt^decade,
                    data = People, cluster = c("decade"))    
results_plot_model <- summary(plot_model)

## get the x01 for the x-axis values in plots
x01 <- min(unique(People$decade), na.rm = T):
  max(unique(People$decade), na.rm = T)

## Branch Effect Over Time
dy.dx1 <- results_plot_model$coefficients["branch"] +
  results_plot_model$coefficients["branch:decade"] * x01 +
  results_plot_model$coefficients["branch:I(decade^2)"] * (x01^2)

dy.dx1_se <- sqrt(results_plot_model$cov.scaled["branch", "branch"] + # V(\beta_{1})
                    results_plot_model$cov.scaled["branch:decade", "branch:decade"] * x01^2 + #x^2*V(\beta_{3})
                    results_plot_model$cov.scaled["branch:I(decade^2)", "branch:I(decade^2)"] * x01^4 + #x^4*V(\beta_{4})
                    results_plot_model$cov.scaled["branch:decade", "branch"] * 2*x01 +
                    results_plot_model$cov.scaled["branch", "branch:I(decade^2)"] * 2*x01^2 +
                    results_plot_model$cov.scaled["branch:decade", "branch:I(decade^2)"] * 2*x01^3)

upr1 <- dy.dx1 + qnorm(.975)*dy.dx1_se # upper bound of the 95 percent
lwr1 <- dy.dx1 - qnorm(.975)*dy.dx1_se # lower bound of the 95 percent

branch_effect <- cbind.data.frame(x01, dy.dx1, lwr1, upr1)
branch_effect$variable <- "branch"

plot_branch <- ggplot(data=branch_effect, aes(x=x01, y=dy.dx1, ymax=upr1, ymin=lwr1)) + 
  geom_line() + geom_ribbon(alpha=0.2) +
  theme_bw() + 
  ylab("Coefficient") + xlab("Year of Birth") +
  ggtitle("Effect of Prominent Branch") +
  geom_hline(yintercept=0, linetype="longdash", colour="red") +
  scale_x_continuous(breaks = c(550, 600, 650, 700, 750, 800, 850)) +
  scale_y_continuous(limits = c(-3.2, 8.5),
                     breaks = c(-2, 0, 2, 4, 6, 8)) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold", size=20),
        axis.title.x = element_text(face = "bold", size=20),
        axis.title.y = element_text(face = "bold", size=20),
        plot.margin=unit(c(0.5, 0.5, 0.5, 0.5),"cm"),
        axis.text.x=element_text(face = "bold", size=18),
        axis.text.y=element_text(face = "bold", size=18)) 
plot_branch


## Exam Effect Over Time
dy.dx1 <- results_plot_model$coefficients["Exam"] +
  results_plot_model$coefficients["decade:Exam"] * x01 +
  results_plot_model$coefficients["I(decade^2):Exam"] * (x01^2)

dy.dx1_se <- sqrt(results_plot_model$cov.scaled["Exam", "Exam"] + # V(\beta_{1})
                    results_plot_model$cov.scaled["decade:Exam", "decade:Exam"] * x01^2 + #x^2*V(\beta_{3})
                    results_plot_model$cov.scaled["I(decade^2):Exam", "I(decade^2):Exam"] * x01^4 + #x^4*V(\beta_{4})
                    results_plot_model$cov.scaled["decade:Exam", "Exam"] * 2*x01 +
                    results_plot_model$cov.scaled["Exam", "I(decade^2):Exam"] * 2*x01^2 +
                    results_plot_model$cov.scaled["decade:Exam", "I(decade^2):Exam"] * 2*x01^3)

upr1 <- dy.dx1 + qnorm(.975)*dy.dx1_se # upper bound of the 95 percent
lwr1 <- dy.dx1 - qnorm(.975)*dy.dx1_se # lower bound of the 95 percent

exam_effect <- cbind.data.frame(x01, dy.dx1, lwr1, upr1)
exam_effect$variable <- "exam"

plot_exam <- ggplot(data=exam_effect, aes(x=x01, y=dy.dx1, ymax=upr1, ymin=lwr1)) + 
  geom_line() + geom_ribbon(alpha=0.2) +
  theme_bw() + 
  ylab("Coefficient") + xlab("Year of Birth") +
  ggtitle("Effect of Keju") +
  geom_hline(yintercept=0, linetype="longdash", colour="red") +
  scale_x_continuous(breaks = c(550, 600, 650, 700, 750, 800, 850)) +
  scale_y_continuous(limits = c(-3.2, 8.5),
                     breaks = c(-2, 0, 2, 4, 6, 8)) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold", size=20),
        axis.title.x = element_text(face = "bold", size=20),
        axis.title.y = element_text(face = "bold", size=20),
        plot.margin=unit(c(0.5, 0.5, 0.5, 0.5),"cm"),
        axis.text.x=element_text(face = "bold", size=18),
        axis.text.y=element_text(face = "bold", size=18)) 
plot_exam



### getting father effects in the same model
dy.dx1 <- results_plot_model$coefficients["FRank"] +
  results_plot_model$coefficients["decade:FRank"] * x01 +
  results_plot_model$coefficients["I(decade^2):FRank"] * (x01^2)

dy.dx1_se <- sqrt(results_plot_model$cov.scaled["FRank", "FRank"] + # V(\beta_{1})
                    results_plot_model$cov.scaled["decade:FRank", "decade:FRank"] * x01^2 + #x^2*V(\beta_{3})
                    results_plot_model$cov.scaled["I(decade^2):FRank", "I(decade^2):FRank"] * x01^4 + #x^4*V(\beta_{4})
                    results_plot_model$cov.scaled["decade:FRank", "FRank"] * 2*x01 +
                    results_plot_model$cov.scaled["FRank", "I(decade^2):FRank"] * 2*x01^2 +
                    results_plot_model$cov.scaled["decade:FRank", "I(decade^2):FRank"] * 2*x01^3)

upr1 <- dy.dx1 + qnorm(.975)*dy.dx1_se # upper bound of the 95 percent
lwr1 <- dy.dx1 - qnorm(.975)*dy.dx1_se # lower bound of the 95 percent

FatherRank_effect <- cbind.data.frame(x01, dy.dx1, lwr1, upr1)
FatherRank_effect$variable <- "FatherRank"

plot_FatherRank <- ggplot(data=FatherRank_effect, aes(x=x01, y=dy.dx1, ymax=upr1, ymin=lwr1)) + 
  geom_line() + geom_ribbon(alpha=0.2) +
  theme_bw() + 
  ylab("Coefficient") + xlab("Year of Birth") +
  ggtitle("Effect of Father's Rank") +
  geom_hline(yintercept=0, linetype="longdash", colour="red") +
  scale_x_continuous(breaks = c(550, 600, 650, 700, 750, 800, 850)) +
  scale_y_continuous(limits = c(-0.2, 0.6),
                     breaks = c(0, 0.2, 0.4)) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold", size=20),
        axis.title.x = element_text(face = "bold", size=20),
        axis.title.y = element_text(face = "bold", size=20),
        plot.margin=unit(c(0.5, 0.5, 0.5, 0.5),"cm"),
        axis.text.x=element_text(face = "bold", size=18),
        axis.text.y=element_text(face = "bold", size=18)) 
plot_FatherRank


grid.arrange(plot_branch, plot_exam, plot_FatherRank, nrow=1)


